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Compact Finite 

Volume Methods for the Diffusion Equation* 


Milton E. Rose 

Department of Mechanical Engineering 
N.C. A&T State University 
Greensboro, NC 27411 


Abstract 

We describe an approach to treating initial-boundary value problems by finite 
volume methods in which the parallel between differential and difference arguments is closely 
maintained. By using intrinsic geometrical properties of the volume elements, we are able to 
describe discrete versions of the div, curl, and grad operators which lead, using summation-by- 
parts techniques, to familiar energy equations as well as the div curl = 0 and curl grad = 0 
identities. For the diffusion equation, these operators describe compact schemes whose 
convergence is assured by the energy equations and which yield both the potential and the flux 
vector with second order accuracy. A simplified potential form is especially useful for obtaining 
numerical results by multigrid and ADI methods. The treatment of general curvilinear coordinates 

l 

is shown to result from a specialization of these general results. 


*This research was sponsored by the National Aeronautics and Space 
Administration under NASA Contract No. NAG- 1-8 12 and by the Air Force Office of Scientific 
Research under Contract No. AFOSR F49620-89-C-0010. 
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Compact Finite 

Volume Schemes for the Diffusion Equation 


Introduction: 

Let V be a domain with boundary surface S on which n is the outward unit normal. 
This paper describes a finite volume scheme for solving the diffusion equation for a potential <$> 
and flu xu in the form 


D.E. 

cjjj = div u - f, 

in V 


\1\ 


grad (}) = u 




B.C. 

a un + = g 

on S, 

t^O 

\2\ 

I.C. 

4> = gQ 

inV, 

t = 0 

\3\ 


in both the steady and unsteady cases. We recall that the solution of this problem satisfies an 
‘energy’ equation 


^ Jv ^ dV + Jv u u dV= Js * uds - Jv * f dV ^ 

which follows by multiplying the first equation in\l\by(J)and integrating by parts. We also recall 
that when f = 0 the maximum- minimum values of(j> either lie on the boundary S or, in case t = 0, 
in V itself. > 

A discrete form of these results will hold for the finite volume scheme and includes 
results about generalized curvilinear coordinate mappings as a special case. We will identify finite 
volume operators div^, grad^, and curl^ which are consistent with their differential counterparts 
and from which a discrete energy expression corresponding to\4\will follow by simple 
summation-by-parts identities. The important identities div^ curl^ = 0 and curl^ grad^ = 0 will 
also remain valid, as we will show at a later point in the paper. A potential form involving the 
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operator div^ gradj^. which is obtained by eliminating u, leads to a symmetric, positive definite 
operator to which multigrid and other fast solution techniques are applicable. 

The fact that the finite volume schemes described here lead to discrete energy 
expressions is the principal result of this paper. It insures that the schemes converge. We will also 
find that the approximations to both 4) and u will be second order accurate. This result is similar to 
that obtained for mixed finite element methods ([1], [8]) and it is possible, in fact, to view the 
present scheme as a finite element method which involves non-conforming elements ([2], [4]). 

Many of these ideas can be illustrated most simply for steady, one-dimensional 
problems, for which reason we first discuss the equations 4> = u, u'= fin detail. We will find it 
natural to introduce a primary mesh, which is formed by the endpoints of subintervals into which 
the basic domain is divided, and a dual mesh, which is formed by the centerpoints of the primary 
mesh. The variables 4> and u will be associated with the primary mesh while another variable 
will be associated with the dual mesh. An algebraic relationship between 4), ip and u on each 
subinterval provides an approximation to the solution operator (for which reason the scheme is 
called compact) and the solution in the large is obtained by imposing continuity conditions for<t> 
and u at points of the primary mesh. Both 4> and qr will be found to converge with second order 
accuracy to the solution at the points at which they are defined. Furthermore, as noted above, 
although u will be defined by one-sided divided differences involving 4> and qr, its convergence 
will also be second order accurate as a result of the continuity conditions imposed. 

In extending these ideas to higher dimensional problems by subdividing the 
solution domain V into volume elements we will also find that a primary and dual grid play a 
natural role. Now the variables 4) and u will be associated with the centerpoints of the faces of a 
volume element while qr will be associated with the center of the element. However, an additional 
variable C (or a box- variable x) which is associated with the edges (vertices) of the element will 
also be required. The compact scheme will describe relationships between these variables which 
produce second order accurate results at the points indicated. In the method discussed here, the 
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non-tangential components of u on a face of an element are obtained with the use of one-sided 
differences in terms of cj) and xy while its tangential components are determined solely by the edge 
values C, and these are obtained in terms of ly at points of the dual grid by bilinear interpolation. In 
the case of (uniform or nonuniform) Cartesian or orthogonal coordinate elements, the variables C 
as well as the tangential components of u on faces of elements can be obtained by postprocessing, 
if required. 

Some of these ideas are familiar from applications of finite volume methods to fluid 
dynamics and are described in Peyret and Taylor [5] and in a review paper by Vinokur [9]. A 
focus of many such methods is on the treatment of conservation laws for inviscid fluids. As noted 
by these authors, the primary and dual grids often play a role in many such schemes whenever 
gradient terms are included, as occur when viscous effects arise. Indeed, the need to accurately 
approximate quantities like the stress tensor at boundaries of general domains is a principal reason 
to resort to finite volume methods, although their use may add significantly to the cost of 
computations. To understand some of the problems which can arise, the diffusion equation can 
serve as a useful model. We will find that the relationships between grid variables differ in 
important respects from those described elsewhere. Many schemes place primary emphasize on 
the vertex (box) variables and most methods deliberately avoid the use of one-sided differencing to 
approximate the flux except, perhaps, at boundaries. Although compact schemes related to those 
described here have been used for Cartesian grids, e.g., [2], [7], the roles of the variables were not 
completely developed. 

The paper concludes by describing the rather staightfoward modifications which are 
required to treat the time-dependent problem. The result is a Crank-Nicholson-type scheme and 
energy arguments provide convergence estimates for the finite volume method. We are also able to 
show that energy arguments can be adapted to a Peaceman-Rachford ADI scheme. 

This paper reflects many valuable insights gained through discussions with 
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colleagues. I mention with particular appreciation T.B. Gatski, D. Gottlieb, H.-O. Kreiss, and 
G. Strang. 
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I. Steady One-Dimensional Problems 


1.1. Notations 

On an interval [1_, 1+] consider the problem : 

D.E. u'= f M.la\ 

= u 

B.C. 4> = g, forx = l ± . M-lb\ 

The associated energy equation is 

J <}>f dx + J" u 3 dx = (4 ju)| + - (<J)u)|_ . M.2\ 

Divide [L, 1+] into M non-overlapping subintervals Ij s {x| xj_ ]/2 ^ x ^ xj+ 1 / 2 } 
with centerpoints xj, j=l/2, 3/2,..., M-l/2. The points xj, j=l,2,...,M-l, are endpoints of the 
subintervals which lie interior to [L, 1+]. We let I e , I c denote the sets of indices corresponding to 
these points: 

I e = { 1 , 2 , . .M- 1 } (interior endpoints) 

I c = { 1/2, 3/2,..., M - 1/2} (centerpoints). 

We adopt the finite-difference notations Axj = xj + j /2 " x i- 1/2» W = Axj u ( x j) = u j. By } 

introducing the central average and difference operators 

M 4>j- ( 4>j+ 1/2 + 4>j_ 1 / 2 V 2 , A 4)j= ( 4>j+ 1/2 - 4>j- 1 / 2 )- 3X 

we can verify the summation-by-parts identities 

A(4>ijj) = (p4)) (Avp) + (A<}>) (mv) M-4a\ 

ji(4>i|») = (ju4>) (pij;) + (A4>) (A»jj)/ 4. M.4b\ 

Both will play a central role in establishing energy results. 
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Figure 1 . An interval showing the association of the variables <{), u with its endpoints and 
of tjj with its centerpoint. 


1.2. A Compact Scheme 

In each subinterval Ij, construct an approximate solution using values 4>j+ 1/2 as 
boundary data. Specifically, 

(i) for j G I c set 


Auj s Uj+ j /2 - uj_ 1/2 = Axj fj M.5\ 

h j u j+ 1/2 = (45+ 1/2 " Vj) 6aN 

h j u j- 1/2 = (Vj - 4>j-l/2) NJ - 6bN 

Adopting the convention that t|r. indicates an approximation to the potential solution at the center of 
the interval, while 4>j±l/2 indicates the approximation at its endpoints we see that M.5\ is a central 
difference approximation to u'= f, while M.6\ approximates 4>'=u by one-sided differences. Then 
eqs.\I.5\, M.6\can be solved for uj+j/2 and ipj in each interval in terms ofchj+i/2 considered as 
boundary data. r 

(ii) Next, require that both u and <J) shall be continuous across every endpoint common to two 
intervals, i.e., 


[u]j = =0 for i € I e (endpoints) M.7\ 

(This is implied by our notation, since we have not distinguished between the right- and left-hand 
limits of u and<J)at endpoints.) Using M.6\ to evaluate Uj + j /2 at the right and left endpoints of the 
adjacent intervals Ii-i/2, Ij+l/2 we 
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M.8\ 


Axj+ 1/2 u i+ 1/2 = (Vi+ 1 - Vi- 1> i € I e 

which indicates that the accuracy of u may be higher than that suggested by the one-sided 
difference expressions which originally defined u inM.6\. Also note that \I.6\ can be used to 
express the continuity condition [u]j = 0 in terms of<J>, ijj with the result 

(<l* " Vi- 1/2 ) / h i- 1 / 2 = (Vi+ 1/2 " 4^) / V i/2 M.9a\ 

and which, in the case of equal intervals, reduces to 

<t>| = pi|;j i G I e (endpoints) \I.9b\ 

The addition of the boundary conditions leads to a determined system of algebraic equations for the 
M values of vy and the M+ 1 values of 4> and u. 

The special nature of this system is worth noting. By eliminating u from M.5\ using 
M.6\we find, in the case of equal intervals, 

(p4>j - ijjj) = 2h^f i e I c (centerpoints) \I. 10a\ 

and if we also write the continuity conditions M.9b\ in the form 

(pijjj - 4>j) = 0 i G I e (endpoints) \L 10b\ 

it is evident that <J)j, i|rj are the odd-even components of the solution vector of a tri-diagonal system 
of equations which is, thus, easily solved. In sharp contrast toM. 10\ the standard finite difference 
scheme for solving M. 1\ is 

(pct>j - i|fj) = h^f i e I c U I e M. 10c\ 

so that some differences with standard numerical treatments of differential equations are to be 
expected using the schemes considered in this paper. 
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1.3. The Energy Equation and Convergence : 


We will now show that this scheme leads to a discrete form of the energy equation 


W\ For simplicity, we will assume a uniform mesh. Add and subtract the expressions 


(Ax/2) uj+i /2 = - (4>j± 1/2 - Vj) 


given inM.6\to obtain 


Ax puj = A4>j 


M.llaN 


vpj = p4>j - Ax Auj/4 


M.l lb\ 


and we see that the approximations to 4> and i|i are second order with respect to the center of an 
interval. Then, starting with 


Auj = Ax fj , M-5\ 

multiply by qx and use M. 1 1\ and \I.4\ to obtain 

Au = (p4> - (Ax/4) Au)-Au M. 12\ 

= A(4>u) - (pu A4> + (Ax/4)(Au)2) 

= A(4>u) - Ax ((pu)2 + (Au)2/4) 

= A((J>u) - Ax p(u^). 

(Here and in the following we omit subscripts when no confusion is likely to arise). Summing 
over centerpoints of intervals 


_ U — 

2.^ ( Vj Auj + Ax puj 2 ) = (<j)u) L ~2.l e 


(M 


Because of the continuity conditions, the jumps [ ] vanish and we obtain the discrete energy 
expression 
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— (U 

2.i c CVjfj+ M u j 2 )Ax = (4 w) l 


M.13\ 


This result will enable us to conclude that <|), 141 and u converge with second order 
accuracy. Introduce centerpoint and endpoint norms as follows: 


I vile — (Xlc Vj ? Ax ) 


1/2 


M. 14a\ 


“ lie = (K +U M + Zle U ?) AX ) 


1/2 


M.14b\ 


By summing M. 8 \ and using M. 6 \ we find 

Vj+i - 4>o = ^0 + 

and if we also assume endpoint conditions (Jjq = 4 )]yj = 0 we obtain the inequality 

II V lie ^ C || U || e . 

The energy expression itself leads to the estimate 

II U HI < II V He || flic 

while the continuity condition M.9\ results in 

II 4> lie — II V He • 

Thus, 

II V lie <C 2 II flic 

II u lie <C II flic 

Il4>l! e <c 2 || f|l c . 


M.15S 


Applying this to the homogeneous problem we conclude that u=({>=v = 0- Thus the 
algebraic problem has a solution which is unique. Next, interpreting f as the truncation error, 
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f = O(Ax^), and we see that<t>, 41 , and u converge with second order accuracy. 


This conclusion about the second order accuracy of u applies not only to its values 
at interior endpoints, as might be expected from M. 8 \ but also to the values uq and u^j at the 
endpoints of [ 1 _, 1 +] itself and which, we recall, were defined in eq. M. 6 \by one-sided difference 


expressions. 

Table 1 presents computations which verify these conclusions for the differential 
equation M.l\ corresponding to the solution <}> = x^(l-x)^. Note that the errors are measured at the 
endpoints of intervals. 



Table 1. 



Error Norms 


# intervals 

II terror lie 

II Uerror lie 

6 

1.234 e-2 

1.851 e-2 

12 

3.472 e-3 

9.259 e-3 

24 

8.680 e-4 

2.893 e-3 

48 

2.170 e-4 

7.957 e-4 

96 

5.425 e-5 

2.080 e-4 

192 

1.356 e-5 

5.312 e-5 


Finally, we remark that more general boundary conditions involving both u and <t> 
also lead to the same conclusions. 


1.4. The Potential Form : 

By eliminating the flux u inM. 1\ we obtain the familiar second order equation (J) = f 
which we call the potential form. The difference scheme which corresponds to this may be 


Compact Finite Volume Methods -(3.2) 
Sat, Sep 16, 1989 


13 


M.E. Rose 


obtained by eliminating u in the scheme just described. To illustrate, using M.6\ eliminate u in 

Auj = Ax fj M.5\ 

to obtain 

h 2 fj = (p<fy - rjTj) , j el c M.16\ 

h = Ax/2 • 

As we have seen, for uniform meshes the continuity conditions [u] = 0 result in 

= fir^. i G I e M.9b\ 

To these are to be added the boundary conditions 

*0 ■ g - '•’m ■ 8+ 

When f = 0 we obtain 

V = H<t> , 4> = 

and the maximum-minimum property of the solution is an immediate consequence, 
used to develop an alternate proof of convergence. 

The potential form for the discrete problem can be seen to arise as a variational 
problem associated with the discrete energy equation. The maximum-minimum property reflects 
the fact that the associated algebraic problem is symmetric and positive definite. In the general 
treatment, the potential form will be more suitable for numerical work, while the general form 
involving 4>, nr and u as a first order system will prove more convenient for theoretical discussions, 1 
particularly for the development of energy estimates. 


M.17N 
This can be 
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II. A Compact Finite-Volume Method 


We now turn to problem of treating the boundary value problem corresponding to 
the steady-state solution of \1\ when V is a general volume. Our objective will be to partition V 
into volume elements 5V in such a manner that the prescribed boundary data on S can be accurately 
transmitted to the boundary elements and then, by solving a discrete boundary value problem in 
each element corresponding to a compact scheme, obtain an approximate solution in all elements 
which also satisfies an energy expression, thus insuring convergence. As a result, we expect to be 
able to treat problems posed in curvilinear coordinates as a special, but important, case. 

In order to be able to generalize the arguments given earlier, a number of additional 
notations will have to be introduced. Before doing so, we will first present a short overview of the 
principal ideas which will be involved. The more detailed discussion and demonstration that the 
final result again leads to an energy expression and thus produces a convergent scheme may be 
omitted if the reader wishes to turn to the discussion of the time-dependent problem given in Part 
III. 

II. 1 An Overview: ) 

Guided by the earlier discussion of the one-dimensional problem, we may identify 
the following requirements to construct a discrete approximation to div u = f, u = grad 4> when 
volume elements are involved: 

(i) construct consistent, discrete approximations to div u, grad <t> in a volume 
element 5V in terms of variables associated with u and 4> at appropriate points of SV , 

(ii) further relate these variables as may be required so as to lead to a determined 
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and consistent system of algebraic equations from which u and <J) can be determined whenever 
boundary conditions are prescribed on S, and 

(iii) develop a difference calculus which allows summation-by -parts identities to be 
applied to volume elements and thereby lead to a discrete energy expression from which a 
convergence argument can be constructed. 


T 



Figure 2. A normal volume element; P, Q,R, are the centers of the volume 5V, a face 5S, and an 
edge on which T is a vertex. Shown at Q is a basis [e] and e ^ , the unit normal to 5S. 


We will confine our attention to hexahedral volume elements 5V, although we will 
allow certain of its vertices to coalesce, thereby including tetrahedra and related elements in our 
discussion. Figure 2 indicates the centerpoint P of 5V, the center Q of an oriented face 5S, and a 
representative centerpoint R of an edge of 5S on which T is a vertex. In addition to the variables 
<KQ), ip(P) introduced in our discussion of the 1 -dimensional problem, we will also associate 
another variable C(R) with points R. These variables will approximate the potential solution at the 
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points indicated. Later, we will also associate a variable x with the vertex points T. Section II.2 
introduces notations which allow us to describe various geometrical properties of an element and 
from which averaging and difference operators, corresponding to\I.4\ can be introduced. 

Section II. 3 describes discrete approximations to the operators div and grad. The 
approximation to div u is based on the use of Gauss’ Theorem in a volume element. Centerpoint 
quadrature approximations to the surface integrals involved leads to a discrete operator divjj which 
provides an intrinsically consistent, second order accurate, approximation to div u in terms of the 
normal components of u at centerpoints of faces 5S. This is a familiar idea in finite volume 
methods. An approximation to u = grad 4> will result from the integral form 

Jc u-dx = 4>(B) - 4>(A) 

in which the points A and B are identified with the points P, Q, and R in Figure 2. A midpoint 
evaluation of the integral will determine two tangential components ofu lying in a surface element 
in terms of the variable £ at edges of5S(Q). The remaining non-tangential component ofu on dS 
will result by using points P and Q and a one-sided evaluation of the integral at Q. This 
approximation will lead to second-order accuracy as a result of imposing additional continuity 
conditions for the variables on 5S, as occurred in the 1 -dimensional case. We indicate by gradh<t> 
the vector u evaluated at each face using this construction. 

The compact scheme which emerges is described in II.4 and can be summarized as 

follows: 

(a) In each element u , 4>, ip, C satisfy divh u = f, gradh4> = u, in which the normal components of 
u arising in div^ u are to be evaluated in terms of the components given by gradh4>. 

(b) Across each oriented face 5S impose the continuity conditions [u-5S] = [<t>] = [Cl = 0. 

(c) If 5S lies on the boundary S of V, then (j) and £ are to be prescribed by the data. 

These conditions exactly mirror the situation discussed for the 1 -dimensional 
problem. However, we shall find that they lead to an underdetermined system of algebraic 
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equations and certain additional conditions will be required. These are furnished by the 
(d) completeness conditions : These conditions determine £ at vertex points in terms of the variable 
ijr by the use of a bilinear interpolation which we indicate as C = ^ (ip). Among other things, This 
will insure that is C bounded in norm by iji. 

This scheme will lead to an energy equation and convergence will be assured. The 
potential form will also show that the maximum-minimum property also holds. 

When curvilinear coordinates are used, many of the geometrical constructions 
required for elements will be provided analytically. Furthermore, for Cartesian or more general 
orthogonal grids, we will find that the variable C can be obtained by postprocessing the results 
using the completeness conditions. This emphasizes the practical advantages of utilizing regular 
elements throughout most of the domain V and restricting the general construction to boundary 
elements. 

II.2. Geometrical Considerations and Notations 
II.2. 1. Normal Volume Elements 

We call {5V} a normal covering ofV by volume elements 5V, themselves called 
normal, if: (i) each element is a hexahedron, some of whose face areas may be allowed to vanish 
providing its volume r remains positive ( degeneracy ), (ii) any face of a boundary element 5V 
which is in contact with the boundary surface S of V is tangent to S at the centerpoint of the face. 
The first assumption has important consequences for the geometrical constructions which will now 
be described. 

Figure 2 indicates a normal volume 5V element with center at P and surface 
elements 5S. Since 8V is a hexahedron, P may be found as the average of the verticies of 5V , 
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while the average of the vertices on a face yields the centerpoint Q and the average of the vertex 
points on an edge yields the centerpoint R. We may thus construct a righthanded, local coordinate 
basis 

[e(P)] = (e 1 (P),e 2 (P) > c 3 (P)) 

at P which is formed by unit vectors in directions which connect the pairs of opposite centerpoints 
of the faces 5S which we denote as Qj+ > i = 1, 2, 3. The face 6S(Qj_) is assumed oriented into 
5V, while 8S(Qj+) is outward; these are sometimes called inflow and outflow faces. We also 
assume that the orientation carries over to neighboring elements by requiring that an outflow face 
from one element be an inflow face on its adjacent neighbor. 

Our primary requirement will be to evaluate the vector u(Q) at the centerpoints of 
faces. In order to do this we will construct a unit basis [e(Q)] at each Q as follows: assume that 
5S(Q) is the outflow face of 5V(P) and the inflow face of a neighboring element 5V(P'), e.g., Q = 
Ql+. We take ej(Q) as a unit vector in the direction of a line connecting P and P'. Also take 
e 2 (Q), e 3 (Q) as unit vectors in directions which connect Q with the centerpoints of the edges R of 
5S(Q) so that [e(Q)]forms a right-handed basis as indicated in the Figure 2. Thus, e 2 (Q), e 3 (Q) 
are tangent vectors on 5S(Q) while e j(Q) is non-tangential. 

We may regard the points indicated as being related by displacement operators Ej+ 
along a direction ej so that 

Qi± = Ej+ P , Ri± j± = Ej+ Qj+ 

with the further understanding that (i, j, k) indicates an even permutation of ( 1 ,2,3). The 
averaging and difference operators defined earlier by \II.3\may be generalized by writing 

Pj <KC) s (<KEi+C) + <KEi_C)V2 , Aj <KC) - (<KEi+C) - <KEi_C)) MI. 1\ 

for a representative centerpoint C. From this it is apparent that the summation-by-parts identities 
M.4\ continue to hold. 

Motivated by this notation we will write also Ajx(Q) to indicate the distance 
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between the points Qj_ and Qj+ and write hj (Q) - Ajx(Q) /2 . Also, h itself will indicate the 
minimum value of hj(P) for all volume elements covering V. 

II.2.2. Covariant and Contravariant Bases on 5V 

A basis [e] = (e j , e 2 , e 2 ) may be considered as a covariant basis from which a 
reciprocal contra variant basis [e] _ l = (e*, eJ, ety can be constmcted satisfying 

e* = (ej* e^) // e, ej = V e(ej* e^) , e'-ej = 5 jj MI.2N 

in which, if ey s ej ej and eU s e^-el, then 

e s det (ey) , 1/e = det (eO). MI.3\ 

Anticipating the connection to curvilinear coordinates, we call ey element metric coefficients at P. 
As a result we may write 

u = £ uje^ = £ u*ei , \II.4\ 

u i = I eliuj , Uj = I eyui 

U U = £ UjU* . 

where Uj = u ej and u' = ue* are the (physical) covariant and contravariant components of u, 
respectively, with respect to [e]. 

Using the unit tangent vectors e2, e 3 the oriented area of a face 5S(Qj +) is 

5S = 1 5S| -e 2 x e 3- = |SS|Vee* MI.5\ 

where 

|6S| = A 2 x(Q) A3x(Q) 

Thus, the covariant component 5S, of 5S is 5Sj = |5S|V e. Also, elementary geometrical arguments 
based on the fact that normal volume elements are hexahedra show that if 

T 0 = A 1 xA 2 xA 3 x 

then the volume measure r of 5V is given by 

T = A iX (Q) ei(PH5S(Qi.) + 5S(Q i+ ))/2 = T 0 MiVe , 
i = 1 ,2,3. If we also assume posiive constants c 0 and cj such that c 0 h^ t ^c t h^ and also set 
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Ve o s PjV e we can summarize these geometrical results as 

T = ip -V eo, 5Sj = Aj xA^ x Ve. 

Finally, we also observe that midpoint quadrature rules result in 
J*udx * u(C) A x(C), 

Ju-dS » u(C)-6S(C) 

Ju dV * u(C) t. 

where C is a representative centerpoint of the figure. 

These facts summarize the principal properties of normal volume elements which 

we shall require. 


MI.6\ 


MI.7N 


II. 3. The Operators div h and grad^ . 

The fundamental geometrical properties of the differential operators div, curl, and 
grad stem from the following: div relates volume integrals with surface integrals, curl relates 
surface integrals with line integrals, and grad relates line integrals with endpoints. We will now 
describe how these properties can be systematically developed for the div and grad operators. 


div h- Gauss’ formula 

J* div u dV = J u-n dS 

when applied on 5V using the quadrature approximations indicated above suggests the definition 

T div^ u = Xj Aj u-5S 

in which the difference operator applies to opposite face centerpoints Qj±. Using the component 
representations for the terms on the right, 

divhu = t -1 £j Aju'-SSj MI 

or, noting MI.6\, 


divhUs(Veo) ^EjAjVeu'/AjX 


MI.10\ 
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grad fr: The approximations to be identified with u = grad <J) will follow from the line- 

integral evaluation 

/ c u-dx = (KB) - <KA) \n. 1 1\ 

in which A, B are endpoints of a curve C. We will describe how the covariant components ofu 
with respect to the bases [e(Q)] can be determined at face centeipoints. 

Tangential Components : 

Recalling earlier notations, suppose e 2 (Qj+), e 3 (Qj+) are tangent vectors to the 
face 5S(Qi+) at Qj+. Evaluating MI. 1 1\ along the line connecting the centers of opposite edges 
R l+, 2 ± and R 1 + 3 + we obtain 

J c u dx = C (Ri+j+) - C (Ri+j-) j = 2,3 

in which we have identified the potential variable on edges of 5 S by C- Next, use the quadrature 
evaluation 


$ C u dx = A j x (u ej(Qi+)) = Aj xuj (Qj+) 
in which uj is the covariant component of u(Q j+). Thus, 

Aj x uj (Qi+) = C (Ri + , j+)-C (Rl+, j-> j = 2,3. 

= Aj C (Qi+) \n.i2\ 

Non-Tangential Components 

We will also require evaluations of the non-tangential components of u at the 

) 

centerpoints of faces and, as for the 1 -dimensional case, want to evaluate these in terms of the 
value of ij/ at the center P of 5V and of <}) at the centerpoints Q of its faces. Let us first introduce the 
one-sided difference operators 


Aj 1 <KQi±) = ± [<KQi±) - v(P)]. 


MI.13N 


Then MI. 8 \ leads to 

I u ds = Aj- 4)(Q i± ) 

in which the line integrals are to be taken between Qj_ and P and between P and Qj+, respectively. 
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We then evaluate the integrals by one-sided quadrature approximations using values of u at the 
faces to obtain 

Judx * (Ajx(P)/2)u(Qi±)-ei(Qj±) 

= (A i x(P)/2)u i (Q i ±) 

so that 

(A i x(P)/2) u i (Q i± ) = Ap (KQi±) = ± l<KQi±) “ V(P)1- ^ 14N 


We can summarize this construction as follows: 

gradh: u(Q) = gradh 4> expresses the set of relationships between the covariant components 

of u(Q) and values of 4>, ip, and C at points Q on the faces of an element given by. 

Tangential components: 

Aj x uj (Q) = Aj C (Q) 

Non-tangential components: 

(Ajx/2) uj (Q) = AptKQ). 

where 

AjitKQ) - ± [<KQ) - ip(P)]- Nn - 13N 


MI. 12\ 


MI.14N 


Eqs.MI. 14\ should be compared to eqs.\I.6V By adding and subtracting these terms 
we obtain an equivalent form: 

'l 

Ajx(P) Pi U}(P) = Aj <KP), 1 5v 

ip(P) = Pj (KP) - Ajx(P) Aj Ui(P) /4 

in which only centered operators at P are involved and which may be seen to be second order 
accurate with respect to P. These forms play an important role in establishing energy results. 
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II.4. A Compact Scheme 

We will now return to the primary problem of solving the boundary value problem 
when V is a general volume with boundary S. Following our earlier discussion, we describe a 
compact scheme for solving the problem as follows: 

(i) in each element 5V solve 

divjj u = f, u = gradjj <J>. MI. 16\ 

Note that divfo u = f provides one equation for the six contravariant components of u on the faces 
of the element, while u = gradj, <J> expresses three equations for the covariant components ofu on 
each of the six sides of 5V, in which the tangential components are evaluated in terms of the 
variable C at the centerpoints of edges while the nontangential components are determined by the 
variables 4> and 4;. 

(ii) use 

u* = £ e’Juj MI.17N 

fromMI.4\to relate the contravariant components to the covariant components arising in MI. 16\.. 

An important simplification will occur in this construction whenever e'J = 0, i*j. 
When this situation arises we call the scheme strongly compact, otherwise we use the term 
weakly compact. For a strongly compact scheme the tangential components ofu on faces will not 
affect the computation directly (and, thus, neither will the variable O but may, if required, be 
determined in terms of the variable C by a postprocessing technique using the completeness ) 

relationships described below. This arises when Cartesian or orthogonal curvilinear grids are 
involved. 

(iii) impose continuity conditions for u. 6, t on faces common to neighboring elements. Incase 
5V is a boundary element, the data on S will be assumed to determine the data 4> and C on any face 
tangent to S by the use of a Taylor series approximation. 
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We are required to calculate u and 4> at the centers of the faces of elements, vjr at the 
centers of the elements themselves, and C on the edges of the faces. By examining a cube, we can 
verify the fact that the number of unknowns involved in the system of equations resulting from the 
steps just described will exceed the number of equations by just the number of interior edges of 
elements, i.e., those edges which are incident to four adjacent elements. This same result will also 
hold for a more general covering {5V} . In order to make the scheme determinate we will add the 
following 

(iv) Completeness Conditions : 

Let N(R) denote the centerpoints P of those volume elements which are incident to 
an edge having R as centerpoint. Using a bilinear interpolation, express C(R) in terms of the 
values ty(R) which are associated with the neighboring elements. The result can be written 

C(R) = ) \n.i8\ 

where the summation includes points P in N(R) and the weights co are non-negative and sum to 1. 
As noted above, when the scheme is strongly compact MI. 1 8\ will allow the tangential components 
of u on faces to be calculated separately. 

II. 5. Ener gy 

We will now show how an energy expression can be derived. From it the existence 
and uniqueness of the solution of the discrete problem will follow. i 

Following the argument given for the one-dimensional case, use the summation 
rules M.4\ and also MI. 1 5\ to obtain 

rdiv^^u = J\Aj(t)u’5Sj 

= £. [jUj4> Aj u' SSi + Aj (J) jjjU 1 5Sj] 

* Lj Kv + 1/4 AjxAjUj) Aj u> 5Sj + Ajxpjuj (pju 1 5Sj)] 

= wZj Aj u 1 5Sj + Ajxpjfuj u* 5Sj) 
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or, recalling that 5S j =/ eAjxA^x, 

r divh 4*1 = <t> r divh u + to £ jMi^i u* V e) 
so that summing over volume elements yields 

ly div^ct*! r= Xy4>divh« T + £y M^u 1 Ve)} to . 

Finally, using the telescoping property of the div^ term on the left-hand side and recalling that 
T = TqV e Q , we find 

Z s 4>u-6S « £y <J> div h u T + £y u u t MI. 1 9\ 

which is clearly consistent with the integral form of the energy equation given by eq.MX. 

Noting that MI. 18\ will insure that C is bounded in norm by i|i, the same arguments as were given 
for the one-dimensional problem will imply that this more general compact scheme converges and 
yields second order accuracy for (b.xu.E and u . 


Compact Finite Volume Methods -(3.2) 
Sat, Sep 16, 1989 


26 


M.E. Rose 



II.6. The Potential Formulation 


MI.20\ 


MI. 2 IN 


We shall indicate how the scheme just described can be formulated solely in terms 
Of*,*, and c, corresponding to the potential form - f . Again following our discussion of 
the one-dimensional case, we write, using MI.9\ and MI.4\ 

T divfo u = Lj Aj u * SSj 

= IiMIj eiju j )5S i 

or, regrouping terms: 

r divh u = IjAj e'Kij 5Sj+ Ij (Aj Xj #i e>juj) 5Sj 

Using MI. 14\ in the first summation we find 

£.Aj e" uj 5Sj = 21^ e" A^ 5Sj/ AjX 

= 4 IjMj e i5 4> 5Sj / Ajx - 4 i|i(P) e U 5Si / Ajx 
in which the summands are evaluated at the centers of faces about the centetpoint P of the element. 
The second group of terms in MI. 20\ involve only tangential components of u at face centets, and 
we may use MI. 1 A to express these in terms of values of C at the edges. In this manner, all of the 
terns inMI.20\may be expressed in terms of the values of*. *, and t in each element. This leads 

to the potential operator div^ grad h 4 >. 

Since 4>, ip, and C are assumed continuous across faces, the continuity conditions 
for u across a face 5Si common to neighboring elements with centers at P and P' reduce to 

[u* 5Sj] = 0. 

The tangential covariant components of u are continuous, since they are determined by C 
Recalling MI. 14\ the condition that the nontangential covariant component 14 be continuous is then 

(AjX + Ajx') ({) = AjX n»'+ Ajx' i|i NIL22N 

which may be compared toM.9a\. 
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Consider the treatment of V^(J) = 0 using a strongly compact scheme. In this case 
e'J = 0, j*i and the condition divjj gradj, 4> = 0 as expressed by using \II.20\ reduces to 

IjMi e 5i (J) 5Sj/ Ajx = n*P) I-iij e ij SSj/AjX MI.23\ 

which may be written 


V = LjOtj cj) \H.24\ 

where 

C4 = (ju i e»5S i /A i x)/(I.p i e i '5S i /A i x), \H.25\ 

Ij«i = 1- 

If we write the continuity condition MI.22\ as 

c{) = Pj i|x + (3j' Pj + Pj' = 1 \n.26\ 

we can conclude that the maximum-minimum property holds. This same conclusion can be 
reached for a weakly compact scheme. 


Referring to the energy expression, it will be seen that the potential form can also be 
viewed as arising from a variational problem. By using the completeness relations, the 
computational problem will thus reduce to determining the variables <j> and i|i associated with a 
symmetric, positive definite matrix and thereby allows a variety of familiar solution techniques to 
be used. 


) 


II.7. Degeneracy 

Recall that the definition of a normal volume element allowed for the area of one or 
more of its faces to vanish. Such degenerate elements play a natural role at the boundary S and are 
especially relevant, as will be shown shortly, in treating problems in curvilinear coodinate at the 
origin. 


Let us suppose that a face with center Q is degenerate, i.e., dS(Q) = 0. Referring 


Compact Finite Volume Methods -(3.2) 
Sat, Sep 16, 1989 


28 


M.E. Rose 



toMI.8\we see that the term involving the factor dS(Q) will not contribute to the evaluation of div^ 
u and, as a result, the contravariant component of u(Q) need not be computed on the face. 

In the potential form, this simply modifies the weight a associated with the degenerate face in 
UI.25\. 

n.8. Curvilinear Coordinates 

We will illustrate how the discussion applies when curvilinear coordinates are used. 
With their use the basis (ej, ej, e^), which was required at each face of a volume element, can be 
constructed analytically, thereby considerably reducing the preliminary computational steps 
required. Using the potential form, we shall also illustrate the treatment of the situation at the 
origin which gives rise to degeneracy. 

We indicate by x = x(q) a mapping in general curvilinear coordinates. The vectors 

dx 

gi = — \n.27\ 

dq 1 

are tangent vectors on the coordinate lines formed by qi = const, and q k = const, so that an element 
of length is determined by 

ds 2 = ££ gjjdq’dqJ 

in which g,j are metric coefficients given by gjj = g, gj. Considering [g] = (gj, gj, gfc) as a 
covariant basis, a reciprocal basis (gi, gj, g k ) is given as 

) 

g 1 = (gj x gk)/^g 

in which g = det(gjj). Also, if gU = g'gj,, then det (gU) = g" 1 . 

Consider a curvilinear volume element dV which is formed as the image of a 
rectilinear volume element whose volume is dq’dqidqK Then the edges of a face on dV given as 
qi = const, are coordinate lines on which gj, gfc are tangent vectors and the oriented area of the face 
is 

dS « (gj x gfc) dqJdq k . 
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We may associate a normal volume element 5V with dV by constructing tangent 
planes through the centers of the surface elements through dV. From geometrical considerations it 
is evident that if we make the following identifications 

Ajxej = gjdqi i = 1,2,3 MI.28\ 

the results of our earlier discussion will apply. If u’(g) indicates a contravariant component of u in 
the curvilinear system we find, in place ofMI. 10 \ 

divjj u ■ (Vgo) “ 1 £. Aj Vg u'(g) / dq 1 . MI.29N 

Similarly, u = grad^ <J>, expressed earlier by MI. 12\ and MI. 14\ now result in 
Tangential components: 

dqj-Uj(g)= AjC MI.30\ 

Non-tangential components: 

(dqV2) uj (g) = AptJ). MI. 3 1\ 


Examples 


The cases of cylindrical and spherical coordinates may help illustrate these results. 
Cylindrical Coordinates : 

Here q = (r, 0, z) and a calculation gives Vg = r and Vg gl 1 = r, Vg g^2 = l/r, 

Vg g33 = r . The centerpoints of faces Qj+ have the coordinates: 

Pl± = (r ± Ai/2, 0, z), ?2± = (r» 0 ± A6/2, z), P 3 + = (r, 0, z ± Az/2) 

} 

Adopting more standard finite difference notations, let <J> (iAr, jA0, kAz) s 4>j j ^ and extend the 

meaning of the spatial averaging and differencing operators p,A to this case by writing (ju , p , 

r 0 


p ), (A , A , A ). We find that in a cell with center at P = (iAr, jA0, kAz) 
z r 0 z 


divfo gradh 4 > = 


4 

r 


{— — t(r + Ai/2) 4^+i/2,j,k + (r - At/2)] 4>i-i/2,j,k 
2Ar 2 

+ -^Me ^ij.k + -^Mz^ij.k - Oq Vij.k} 
rA0 Az 


MI.32N 
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in which 


I 


o 0 =r( 


1 


1 


Ar 2 Az 2 


■) + ■ 


rA0 


Assuming equal spacings, the continuity conditions across faces have the form 

4 > = MV 


\II.33\ 


in which p = (ju , u , u ). 

r 0 z 

A cell with center at (Ai/2, jA0, kAz) has the origin as a degenerate face. In this 
case we see that the coefficient of the value of the unknown <J}q j ^ vanishes, so that no assumption 
other than boundedness at the origin is required for the scheme. 

For cells not involving boundary points, we may use the continuity condition in 
MI.29\to express the result solely in terms of ip. By adopting the notation s A f / Ar , etc. to 

indicate divided differences, the potential operator for such cells reduces to the form 


div^ gradh 4> = { &r + + ^ $e + ^ z } ^ MI.34N 


which can be seen to be consistent with the familiar differential form in cylindrical coordinates. 
However, we emphasize again that the boundary conditions involving <J> cannot be handled directly 
from this form but must, instead, be developed using MI. 30\ directly. 


) 


Spherical Coordinates: 

In this case q = (r, cp, 0) and a calculation shows Vg = r^ sin cp while 
VggH = r 2 , Vgg 22 = sin cp , Vg g 22 = 1. Again assuming a uniform mesh spacing, we find in 
an element with center at P = (r, cp, 0) that 
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divh gradfo <t> = — — — { — [(r + At/2) <t^+i/ 2 .j.k ( r — Ar/2) 4>i-i/2j,k 1 

n r 2 sin X l 2Ar 2 

+ — - — [sin(x + Ax^^y+i^.k + sin(x - Ax^M^j-i^.k! 
2rA X 2 

e^i.j.k " °0 Vi,j,k } 

A0 

\C.35\ 


in which 

Oo = { — (r 2 + Ar 2 /2)sinx + — L- sinx cosAx/2) + } • 

Ar 2 Ax 2 A0 


With |i = (ju , fj , ji Q ) the continuity conditions again lead to <t> — jaiy. We see that 
r X 0 

the situation at the origin is exactly similar to that in the previous example; in particular, only the 
boundedness of the solution at the origin is required. Away from boundary cells, the continuity 
conditions can be used to reduce matters to the form 


A t ip = — ^ — { (r 2 + Ar 2 /2) sinx 5f + 2 sinx 5 r Mr 
r^inx 

+ sinx cos A X /2)5^ + 5e } V 


\D.36\ 


which can be seen to be consistent with the potential form of the diffusion equation in spherical 
coordinates. We again add a precautionary note about attempting to treat the boundary conditions 
involving 4> from this form. 

II.9. The curi^ Operator 

In our treatment of u = grad^ cf), emphasis was placed on evaluating the tangential 
components of u at the center Q of each face 5S using a covariant basis [e(Q)]. Using MI. 1 A 
these were expressed in terms of the variable C which was associated with the centerpoints of the 
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edges of the face. The completeness conditions MI. 18\ were then used to relate C to iy. The 
contravariant component of u which arose in the definition of divjjU were evaluated in terms of 
these covariant components using MI.4\. 

It is desirable, also, to also have a finite volume approximation to the curl operator 
which we will indicate by curl^u. A satisfactory calculus would then yield both div^ curl^u = 0 
as well as curl^ grad^ ({> = 0 as identities. We shall now indicate how this may be accomplished. 

In Figure 2, T indicated a representative vertex of an element. Let [e(T)] indicate a 
basis formed by the edges which intersect at T using the orientation established by [e(P)]; also let 
e(R) indicate the translation of e(T) to the centerpoint R of an edge TT\ A suitable definition of 
curljj is suggested by the following: If w = curl u, then a quadrature approximation in Stokes’ 
formula 


leads to 



urn dS = u-ds 


00-5S “ x-Aj (u-efc) - Aj x-A^ (u-ej) 

involving the edges of a face 5S. Using component representations, we are led to define curl^ in 
terms of th e set of contravariant components to 1 associated with the faces of 5V by means of 

to = curlh u =*• co 1 s (Aj ufc /Aj x - A^ uj /A^ x) MI.32N 

in which the covariant components u = u -e. With this definition the identity 

divjj curljj u = 0 MI.33N 

follows by noting the cancellation of ‘line integrals’ along the common edge of faces. (It is also 
possible to identify a vector associated with the centerpoint of dV, say <curlhu> s LPito‘ej(P), 
but this vector will not be annihilated by divj, rind for which reason will not be considered further 
here). 


As already noted, our earlier treatment of u = grad^ (j> will not furnish a direct 
means of evaluating the covariant components u occurring above and thus we cannot interpret curljj 
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grad h <t>. In order to overcome this, introduce a new variable x which will be defined at vertices T 
and is often called a box-variable. Recalling the intrinsic definition of grad indicated by \II. 1 1\ 
we will define gradjj along edges by 

A i xu i (R) = A i x(R) VI1 - 34 

in which Ajx is the distance between vertices along an edge with direction ej on which R is the 
centerpoint andAjX (R) is the difference in values of x at the vertex neighbors of R. Then 

U) i (a) S (A j A k x -A k A jX ) * Aj x-A k x s 0 


i.e. 


t 


curl k grad k 4> s 0. 


MI.35\ 


This suggests a theoretical advantage in using the box- variable x instead of C which 
was used in the earlier discussions. The most direct way to accomplish this is to define the edge 
variable C(R) as the average of x at vertices neighboring R, i.e., C = MX. and also re-express the 
consistency conditions MI. 1 8\ so that now x(T), rather than C(R), will be related to vjt(P) at each 
interior vertex by a bilinear interpolation of the form x = H(ijf), which will now involve the six 
neighboring elements having T as a vertex. 
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III. The Time-Dependent Problem. 


III. 1 The Potential Form . 

The potential form of the time-dependent diffusion equation \1\ is, with f = 0, 

4^ = div grad $ Mil. 1\ 

The discrete finite volume operator divjj grad^ <}> corresponding to div grad 4> has been shown to 
yield a discrete energy expression which corresponds to the time- independent energy terms in\4\. 
The principal question remaining, therefore, is how to include the term in the numerical scheme 
which will provide the approximation to the term J <J> 2 < dV in the energy expression. 

Divide [0,T] into N equal parts of width At, write ^ = n At , and adopt the standard 
finite difference notation u(nAt) = u n . Also, let p t , Aj indicate the central averaging and difference 
operators which effect the time index as indicated: 


jj j u n = (u n+ ^ + u n fl 

Aj u n = ( U n + l/2 - u n-l/2^ 


MU.2\ 


From M.4\ it is evident that 


p t vji n A t i|t n = T A t (i|t n )2. 
We will consider the following three-level 
A. Compact Scheme . 


MII.3\ 


Aj vjj n = At divjj gradjj cj/ 1 MII.4a\ 

ip n = MII.4b\ 

to which the appropriate continuity and completeness conditions MI. 22\ and MI. 18\ are understood 
to apply. By multiplying MII.4\by vp n we obtain 

1 A t (n; n )2 = At Hf n divfo grad^ cj) 0 MII.6\ 

2 

in a volume element . The term on the right-hand-side has already been shown to lead to the 
energy expression MI. 19\. Thus, by summing MII.6\ over elements in a time-strip and using the 
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implied continuity conditions across faces we will obtain 

Xy 1 A t (ijj n ) 2 T + At Xy u n «m t = At £ s <J> n u n -5S \m.7\ 

which is the discrete energy expression which was expected. 

Assuming dissipative boundary conditions in the sense of Kreiss [3], it is evident 
from this result that we can conclude: both u/ n and 6 ° are bounded in a discrete norm by the initial 
and boundary data. Because the scheme is consistent, this also implies that the scheme converges. 
(Our previous discussion of the maximum-minimum property can allow us to conclude — 
convergence as well.) 


We can summarize this in several equivalent ways, in each of which the additional 
continuity and completeness conditions are assumed to apply: 


B: 2-step : 


V n+ 1/2 _ = (At/2) divfo gradfo 4> n+ 1 / 2 \III.8a\ 


and 


V n+ 1 - y n+ 1/2 = (At/2) divjj gradjj <J>n+ 1 /2 \III.8b\ 

As above, this is a two-step scheme in which the first requires solving the equations in a time strip 
for <j> n+ 1/2 and l | f n+ 1/2 implicitly in terms of the boundary data for <Kt n +l/2) and the initial data 
y n . With these values determined, then ty n+ l can be determined explicitly from the second set of ,> 
equations and furnishes initial data for the next time strip. 

C: 1-step: 

Eliminating ly 11 in the above scheme we find 

A( y n = At p t divh grad^ 4>n MII.9N 

This is a one-step Crank-Nicholson form . 
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Example : 


Using our earlier 1 -dimensional notations, consider the treatment of 4^ = 4>xx . 

We find that the 1-step scheme may be written in potential form as 

Aj ij/j n = k p t (p4>j n - i € I c (centerpoints) ME. 10a\ 

where k = At /h^, h = Ax/2 To this we add the continuity conditions in the form 

jjHrj n+ l/2 _ (jj.n+1/2 = o i g I e (endpoints). MH.lObN 
Comparing withM.9\we again recognize that the variables 4^ t|fj are the odd-even components of a 
tri-diagonal system and are therefore readily obtained at each time step. 

III.2 An ADI Scheme : 

Except for the one-dimensional example just considered, the basic scheme Mil. 4\ is 
implicit and effective solution methods must be considered in order to treat the general problem in 
higher dimensions in a practical manner. It is especially desireable that a proposed method also 
provide an effective means of treating the steady-state problem as well. We will now show that a 
Peaceman-Rachford-type ADI scheme, using a straightfoward treatment of intermediate boundary 
conditions, can solve the finite volume problem with second order accuracy in both space and time. 
Again, the existence of an energy estimate will provide the key to convergence. 

We will find it convenient to make a few slight changes in some of our notations. 

First, we will let v indicate the time average of a variable v, i.e., let 

v 11 = p t v n MIU1N 

Also, we let 

Fj n = T^AjCuH-dSj) MII.12N 

so that MI.9\ leads to 
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div h u n s F 1 n + F 2 n + F 3 n . 

Then, for example, the 1-step scheme MIL ION can be written as 

A t iy n = At (F i n + F2 n + F3 n ) 
in which u = grad^ (j>. Finally, we indicate by g n assumed boundary data for <{) on S at time t n . 

Limiting our discussion to two-dimensions, consider the following 2-step scheme: 
v n+ 1/2 _ = (At/2) (F ! n+ 1/2 + F 2 n ) Mil. 14a\ 

t|» n+ l - ij, n+ l/2 = (At/2) (Fi n+1/2 + F 2 n+1 ) Mn.l4b\ 

and note its similarity to MII.8\ Recalling our discussion of the one-dimensional example, it is 
evident that the first step cam be solved as a one-dimensional problem using boundary conditions 
4» n+ ^ 2 = g n+ ^ 2 . We solve the second step similarly, but with the boundary conditions 
4> n+1 =g n+1 . . 

Adding and subtracting MIL 14a\andMII. 14a\ and relabeling the time level for 
purposes of later discussion, we find 

v n+ 1/2 _ v n-l/2 = At (Fj n + F 2 n ) MIL 15a\ 

V n - y n = (ai/ 4) ( F 2 n " m - F 2 n+ 1/2 ). Mn. 1 5b\ 

The consistency of (a) with the 1-step Crank-Nicholson form MIL 1 1\ is evident, the truncation 
error being 0(At 2 X Eq.(MII. 15b) is also consistent with MIL 4(b)\ since the right hand term is 
0(At^X We may regard MIL 1 4\ either as a scheme to approximate the solution of the compact 
scheme MII.4\ or as an independent scheme to solve the differential equation. / 

We shall show that this scheme leads to an energy estimate which is consistent with 
MILA to within terms of 0(At 2 ). Recalling the basic energy argument inM. 1 l\and the definition 
of Fi in MIL 1 2\ we find 

ij/ Fl = A((Jxx^dSj) - pjCu^i) 

V • F 2 = A(4>u 2 dS 2 ) - p 2 (u 2 u 2 ) 

so that by multiplying M II. 1 5a\ by ly and then using MIL 1 5b\ we find 
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— AfOy 11 ) 2 t + [pj(u^uj) n +p 2 ( u 2 u 2 ) n ] = [a(<JxjMSj) + A(<Jxi 2 dS2)] + 0(At 2 ). 

2At 

Summing over volume elements and noting that ijj = 41 + Q(At 2 ), etc., we again obtain the energy 
expression MIL A to within terms 0(At 2 ). This result is sufficient to enable us to conclude that the 
ADI scheme converges and furnishes an approximation with the same degree of accuracy as the 
compact scheme MH.4N. This result is independent of any of the fixed ratios At / (Axj) 2 occuring in 
either scheme. 

To treat the 3D problem, consider in place ofMII. 14a\ etc., the following: 

v n+l/4 _ v n = (At/4) (2Fi n + F2 n+1/4 ) MILI6S 

i|/ n+ l/2 _ ^n+1/4 = (At/4) (F2 n+1/4 + 2F3 n+ly,2 ) 

^n+3/4 _ ^n+1/2 = (At/4) (F 2 n+3/4 + 2F3 n+1/2 ) 
v n+l _ ^n+3/4 = ( At /4) (2F 1 n+1 + F 2 n+3/4 ). 
to which the intermediate boundary conditions 4> n+ ^ = g n+ *S k = 1/4, 1/2, 3/4, 1 apply. Then 
ipn+1 - v n = At [(F 1 n+1 + Y { n )/2 + (F 2 n+3/4 + F 2 n+1/4 )/2 + F 3 n+1/2 ] 
so that the scheme is consistent with the Crank-Nicholson form. Also, 

( v n+l + ¥ ny 2 _ ^n+1/2 = ( At /4) [(F 1 n+1 - Fj n ) + (F2 n+3/4 - F2 n+1/4 )]. 

The last pair of equations, which correspond to MIL 15\in the 2D case, can provide the basis for 
an energy argument, although we omit the details here. 


} 
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Conclusions: 


We have described a finite volume method which closely maintains the parallel 
between differential and difference arguments. By using intrinsic geometrical properties of the 
elements, we were able to describe discrete versions of the div, curl, and grad operators which led, 
using formal summation by parts techniques, to discrete energy equations as well as to the 
identities divjj curl^ u = 0 and curl^ gradfo 4> = 0. The solution of the initial- boundary value 
problem for the diffusion equation was described directly in terms of these operators by compact 
schemes and the resulting energy equations insured convergence. The schemes could also be 
simplified to a potential form which can offer computational advantages. Finally, the treatment of 
general curvilinear coordinates was shown to result from a specialization of these results. 
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